The present application claims priority on European Patent Application 00306780.8, filed on Aug. 9, 2000.
The present invention relates to a method of processing an image, for example a seismic image, applying a diffusion process. Such methods are known, see for example the thesis J. Weickert, xe2x80x98Anisotropic Diffusion in Image Processingxe2x80x99, Ph.D. Thesis, University of Kaiserslautem, Germany, January 1996, pages 41-44.
Known is a diffusion method that comprises the steps of
(a) obtaining an n-dimensional initial image data set, wherein each element u(m=0) of the data set is the initial image intensity of a point of the image, and wherein n=2 or n=3;
(b) calculating for each point the partial derivatives of u(m) in n directions to obtain n derivative data sets of the partial derivatives uxi;
(c) calculating for each point a symmetric nxc3x97n structural matrix S, wherein the elements spq equal uxp.uxq;
(d) calculating for each point an nxc3x97n diffusion matrix D, wherein the elements dij are a function of the elements spq;
(e) calculating for each point u(m+1) from u(m) using the following equation u(m+1)=u(m)+div(Dgrad(u(m))); and
(f) repeating steps (b) through (e) M times to obtain the processed image.
Unless otherwise specified, the following holds throughout the specification and the claims:
i, j, p, q=1, . . . , n ;
uxi=∂u/∂xi.
An element of the image data set is also referred to as a point of the image.
In the known method an initial image is processed to obtain successive images, wherein each image is obtained from the previous one by a diffusion process. The diffusion method can as well be described by the       diffusion    ⁢          xe2x80x83        ⁢    equation    ⁢          xe2x80x83        ⁢                  ∂        u                    ∂        t              =            div      ⁡              (                  D          ·                      ∇                          xe2x80x83                        ⁢            u                          )              .  
There are four types of diffusion methods. At first there is the isotropic diffusion in which the diffusivity D is a scalar and the anisotropic diffusion in which D is a matrix. Then for each of these two types, the diffusion can be linear, in which D is not a function of u or the diffusion can be non linear in which D is a function of u.
This method can as well be applied to processing a seismic image, in that case the image intensity is the amplitude of the seismic signal.
An advantage of the diffusion method is that noise is suppressed. A further advantage, which is particularly relevant to processing a seismic image is that after a few cycles the geological structure is highlighted more clearly. However, a problem with this method is preserving edges.
Edge preservation is of great interest in the interpretation of seismic images, which play an important role in the study of underground formations, and in particular in the study of underground formations that contain hydrocarbons. Moreover there is a great interest in edge preserving techniques, which allow highlighting of faults while removing noise from the seismic data.
For edge preservation in combination with isotropic diffusion processing is proposed a method based on the Perona-Malik model (see the thesis of Weickert, page 10) of which the diffusion equation is                     ∂        u                    ∂        t              =          div      ⁡              (                              g            ⁡                          (                                                "LeftBracketingBar"                                      ∇                                          xe2x80x83                                        ⁢                    u                                    "RightBracketingBar"                                2                            )                                ⁢                      ∇                          xe2x80x83                        ⁢            u                          )              ,
wherein the non-linear diffusivity is g(|∇u|2)=(1+|∇u|2/xcex2)xe2x88x921, the scalar xcex being a pre-determined edge preservation parameter.
Applicant proposes a simple addition to the known method that has proved to be an effective way of preserving edges while capable of removing noise. Moreover Applicant proposes such a method that is applicable to anisotropic diffusion.
To this end the method of processing an image according to the present invention comprises the steps of
(a) obtaining an n-dimensional initial image data set, wherein each element u(m=0) of the data set is the initial image intensity of a point of the image, and wherein n=2 or n=3;
(b) calculating for each point the partial derivatives of u(m) in n directions to obtain n derivative data sets of the partial derivatives uxi;
(c) calculating for each point a symmetric nxc3x97n structural matrix S, wherein the elements spq equal uxp.uxq;
(d) calculating for each point an nxc3x97n diffusion matrix D, wherein the elements dij are a function of the elements spq;
(e) calculating for each point u(m+1) from u(m) using the following equation u(m+1)=u(m)+c.div(xcex5.Dgrad(u(m))), wherein c is a predetermined constant, 0xe2x89xa6cxe2x89xa61, and wherein xcex5 is a scalar, 0xe2x89xa6xcex5xe2x89xa61, wherein xcex5 is close to zero when near an edge and close to 1 when far away from an edge; and
(f) repeating steps (b) through (e) M times to obtain the processed image.
In this way edge preservation is introduced in anisotropic, non-linear diffusion.
Suitably step (c) comprises calculating for each point a symmetric nxc3x97n structural precursor matrix S0, wherein the elements s0pq equal uxp.uxq; creating (xc2xd)n(n+1) data sets, wherein the elements of the first data set are s011 pertaining to each point, the elements of the second data set are s012 and so on; and filtering each of these data sets by convolution with a suitable kernel to obtain the elements spq of the nxc3x97n structural matrix S.
Alternatively step (c) comprises filtering for each point the partial derivatives uxi by convolution with a suitable kernel to obtain regularized partial derivatives urxi; calculating for each point a symmetric nxc3x97n structural precursor matrix S0, wherein the elements s0pq equal urxp.urxq; creating (xc2xd)n(n+1) data sets, wherein the elements of the first data set are s011 pertaining to each point, the elements of the second data set are s012 and so on; and filtering each of these data sets by convolution with a suitable kernel to obtain the elements spq of the nxc3x97n structural matrix S.
The above described alternatives for step (c), wherein a precursor matrix is calculated allows defining the edge preservation parameter as follows: xcex5=trace(S0S)/(trace(S0).trace(S)). In this way the edge preservation parameter is calculated for each point.
The present invention also relates to a method of processing an image in order to display the edge preservation parameter. This method comprises the steps of
(a) obtaining an n-dimensional initial image data set, wherein each element u(m0) of the data set is the initial image intensity of a point of the image, and wherein n=2 or n=3;
(b) calculating for each point the partial derivatives of u(m) in n directions to obtain n derivative data sets of the partial derivatives uxi;
(c) calculating for each point a symmetric nxc3x97n structural matrix S0, wherein the elements s0pq equal uxp.uxq; creating (xc2xd)n(n+1) data sets, wherein the elements of the first data set are s011 pertaining to each point, the elements of the second data set are s012 and so on; and filtering each of these data sets by convolution with a suitable kernel to obtain the elements spq of the nxc3x97n structural matrix S; and
(d) calculating for each point xcex5=trace(S0S)/(trace(S0).trace(S)) and attributing the calculated value of xcex5 to each point to obtain an image in which faults are highlighted.
Alternatively, the method of processing an image comprises the steps of
(a) obtaining an n-dimensional initial image data set, wherein each element u(m=0) of the data set is the initial image intensity of a point of the image, and wherein n=2 or n=3;
(b) calculating for each point the partial derivatives of u(m) in n directions to obtain n derivative data sets of the partial derivatives uxi;
(c) filtering for each point the partial derivatives uxi by convolution with a suitable kernel to obtain regularized partial derivatives urxi; calculating for each point a symmetric nxc3x97n structural precursor matrix S0, wherein the elements s0pq equal urxp.urxq; creating (xc2xd)n(n+1) data sets, wherein the elements of the first data set are s011 pertaining to each point, the elements of the second data set are s012 and so on; and filtering each of these data sets by convolution with a suitable kernel to obtain the elements spq of the nxc3x97n structural matrix S;
(d) calculating for each point xcex5=trace(S0S)/(trace(S0).trace(S)) and attributing the calculated value of xcex5 to each point to obtain an image in which faults are highlighted.
In the paper xe2x80x98Coherence-Enhancing Diffusion Filteringxe2x80x99 by J Weickert, International Journal of Computer Vision vol 31(2/3), pages 111-127, (1999), modifications of the method known from the abovementioned thesis are discussed. The paper discloses a method for enhancing coherence in images with flow-like structures. This method is not a method for preserving edges, and no parameter corresponding to the parameter xcex5 of the present invention is disclosed.